The effect of quantum fluctuations on the coloring of random graphs 
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We present a study of the coloring problem (antiferromagnetic Potts model) of random regular 
graphs, submitted to quantum fluctuations induced by a transverse field, using the quantum cavity 
method and quantum Monte-Carlo simulations. We determine the order of the quantum phase 
transition encountered at low temperature as a function of the transverse field and discuss the 
structure of the quantum spin glass phase. In particular, we conclude that the quantum adiabatic 
algorithm would fail to solve efficiently typical instances of these problems because of avoided level 
crossings within the quantum spin glass phase, caused by a competition between energetic and 
entropic effects. 
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^ , Among the fascinating potential applications of quantum mechanics, the promise of quantum computing is to make 

' use of the laws of quantum mechanics to enhance computers' calculation power. Besides the great effort of research 
towards the physical realization of quantum computers, a lot of activity has been devoted to develop quantum 
algorithms and understand their efficiency. A central problem encountered in almost all branches of science is to 
optimize irregularly shaped cost functions: the Quantum Adiabatic Algorithm (QAA) is a generic and universal 
• procedure to tackle such problems. Suppose one wishes to find the ground state of a Hamiltonian Hp acting on N 
^ ' qubits. The idea of the QAA is to implement an interpolating Hamiltonian H(r) ~ Hp+THq such that the quantum 
I ' . ' computer can easily be initialized in -ffq's ground state, and to decrease the interpolation parameter F from a very 
Ci , large value down to zero. If this is done slowly enough the adiabatic theorem ^] ensures that with high probability 
(/3 ■ the system remains at all times in the ground state of the interpolating Hamiltonian. In particular, at the end of 
the evolution, it is in the ground state of Hp and the original problem is solved. The crucial question is how slow 
the evolution should be in the thermodynamic limit iV — > oo, and whether the total evolution time T has to grow 
polynomially or exponentially fast with N. This is reminiscent of the classical complexity classes 0, which indeed 
have quantum counterparts |9l-lll| . Quite generally, the adiabaticity condition requires T to be larger than the inverse 
^ , of the squared gap between the ground state and the first excited state of H(r), for all F. Hence, the efficiency of 
Q ' the QAA mainly depends on the rate of closing of the minimal gap along the interpolation, even though more subtle 
issues such as determining the residual energy after the interpolation [l2l [isj would require a further understanding 
of the annealing dynamics. 

^-H ' The possibility that the QAA could outperform classical algorithms has triggered a lot of work on optimization 
^ , problems in quantum fields (see [l3 - [T7| for recent reviews), trying in particular to pinpoint the rate of closing of 
' the minimal gap along the interpolation, path. As it has now become usual, these studies used random ensembles 
of constraint satisfaction problems [isl . [l9l |. focusing on results valid with high probability in the thermodynamic 
^5 limit. Early results gene rated excitement by reporting polynomial scaling of the minimal gap for some classically 
' hard problems [sl. [20n2]|: however, these results where hampered by strong finite-size effects or by the fact that the 
instances considered were not typical of the underlying problem. On the other hand, [23-125'] obtained negative results 



for the success of the QAA on a certain class of optimization problems: namely, they found a first order quantum phase 
transition during the interpolation between the quantum and the classical Hamiltonians. Such a phase transition is 
" . . ' not surprising in the context of fully-connected quantum spin gl asses (26l - l29| , and is known to generically lead to a 
^ , gap vanishing exponentially fast with the system size N j23 - l24 | and thus to a blow-up of the time needed by the 
k>( QAA. Failure of the QAA because of another mechanism was also discovered for different models in [30l - l32l |: because 
5h ' classical energies have a non-trivial perturbative expansion in F, the energy of an excited state may decrease 

d ■ much faster than the one of the ground state, leading to avoided "perturbative" crossings near the classical end of the 
algorithm. However, these results were obtained either on a toy model, or using well-chosen Hamiltonians Hp with a 
peculiar structure of the low-energy spectrum of i?p, that were not typical of the optimization problem considered. 

Henceforth, one would like to understand what happens more generically for models with multiple and not nec- 
essarily isolated ground states. This is expected to happen for typical optimization problems, which are known to 
possess complex and intricate configuration spaces that were unveiled by the use of statistical mechanics tools ^33l - f36j . 
In these more common cases, do the avoided crossings remain finite and isolated in F (hence leading to singularities in 
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the ground state energy for N — oo) or do they proHferate and accumulate, leading to a continuum of level crossings 
and a gapless phase? Hints for an answer have been obtained on a toy model in [itI . [3^ ; this paper aims to highlight 
their conclusions on a more realistic optimization problem. This question is also important from a practical point 
of view as it has been argued that a finite number of level crossings could be eliminated by a suitable redefinition 
of the quantum Hamiltonian j37l - l40| . Also note that the Xorsat problem studied in [24] already possessed in some 
cases an exponential degeneracy of its ground states, but they had a particular structure that smoothened the effect 
of quantum fluctuations. 

The model that we shall study is the coloring one, a famous problem in combinatorics, which is known to be 
classically very hard to solve - more precisely, NP-complete [4l|. Given a graph and gcoi colors, it consists in coloring 
the vertices such that no edge connects two vertices of the same color. In physical terms, this corresponds to a 
(7coi-states antiferromagnetic Potts model. This model can be studied on any kind of graph, but we shall focus in the 
following on random regular graphs, with four colors (gcoi — 4) and connectivity c — 9 and c — 13. Thanks to the 
quantum cavity method [4^ l43j , we will compute the phase diagram of this model in both cases and unveil the nature 
of its quantum spin glass phase. Our main results are (i) the nature of the quantum phase transition that occurs 
along the interpolation from the quantum phase to the classical one at low temperature, that we find to be continuous 
in the spin glass language and thus of third order thermodynamically, and (ii) the presence of a continuum of level 
crossings within the spin glass phase, that are induced by entropic effects and in particular by the clustered structure 
of the spin glass phase. This last feature should lead the QAA to fail to solve the coloring problem efficiently. Our 
results will also be corroborated by Monte-Carlo simulations. Note that result (ii) confirms the predictions made in 
(321 1 on a much simpler toy model. 

The plan of the paper is as follows. We first briefiy recall in Section the classical features of the model and the 
definition of its quantum version. Section IIIII presents the phase diagrams that we obtain for the quantum coloring 
problems, while Section ITVl presents in greater details the structure of its quantum spin glass phase. Technical details 
are deferred to a series of Appendices. A brief part of our results (mostly Sec. IIV A[) already appeared on the review 
paper (l7l |. 



II. THE COLORING PROBLEM: CLASSICAL PICTURE AND DEFINITION OF ITS QUANTUM 

VERSION 



A. The classical coloring problem 

Let us consider a graph G = {V, L) with V a set of N vertices and L a set of Af edges between pairs of vertices. 
We introduce a Potts variable di e {1, . . . , gcoi} on each vertex i of the graph, and denote g_— {ai, . . . ,(Jm) the global 
configuration of the variables. To each of these configurations we associate an energy, or cost function, 

E{<l) = 5,.,.^. , (1) 

where the sum runs over the edges of the graph, and 5 denotes the Kronecker symbol. Interpreting the Potts variable 
(Ji as a color given to vertex i, this cost function counts the number of monochromatic edges in the configuration g_. 
In physical terms there are antiferromagnetic interactions between pairs of variables linked by an edge of the graph. 
The Gibbs-Boltzmann probability measure at inverse temperature /3 is then defined as 

with the partition function Z{[3) ~ ^-P^is.) ensuring its normalization. 

As mentioned in the introduction the coloring problem is NP-complete (for q^oi > 3) in its decision version, namely 
there does not exist any algorithm able to decide the existence of a proper coloring (a configuration of zero energy) 
in a number of operations bounded by a polynomial in M for every graph G. Computer scientists ^19l] thus turned 
to random ensembles of graphs in an attempt to study the typical difficulty of the coloring problem, typical meaning 
here "with a probability going to one in the large size (thermodynamic) limit TV — )■ 00" . The interesting regime is the 
one of sparse random graphs, with the number of edges M of the same order than N . The absence of an underlying 
Euclidean structure in these random graphs makes the correspondin g st atistical model of the mean-field type, hence 
the methods devised for the study of mean-field spin glasses models [4J could be applied to characterize the typical 
behavior of t he p artition function Z and of the measure /x. For the coloring problem, see [i^ for the use of the replica 
method and [35|, I46l - l49j for the application of the cavity method. The main outcome of these studies is the unveiling 
of phase transitions in the behavior of the free energy density f — — In Z (that concentrates around its average in 
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the thermodynamic hmit) and in the properties of /x, as a function of the inverse temperature p and of the parameters 
of the random graph ensemble. 

In the following we shall sketch these phase transitions and the cavity methodology used to derive them; for more 
details and justifications of the equations the reader is refered to the above quoted references for the coloring problem, 
and to [sgI. l5Cll| for more generic presentations of the cavity method. We will focus for technical simplicity on the 
case of random c-regular graphs: the graphs G are chosen uniformly at random among all the graphs on N vertices 
in which all vertices have the same degree (number of neighbors, or connectivity) c. In the thermodynamic limit 
(A'^ — >■ oo), such graphs are locally tree-like, meaning that if one selects one of its vertex at random, the shortest loop 
around it is larger than any fixed length with a probability which goes to one when N ^ oo [sij . The cavity method 
exploits this tree-like character, building on the exact solution for finite trees that is easily obtained by recursion, 
and dealing with the boundary conditions induced by the long loops of the random graphs in a self-consistent way. 
Consider an arbitrary variable i in a large random c-regular graph, and its marginal probability r]{(7i) that would be 
obtained from the Gibbs-Boltzmann measure ^ if one of the c edges around it were removed. Forgetting the possible 
effects of the long loops of the graphs, the translational invariance and the recursive structure of the tree implies that 
77 is solution of the following self-consistent equation, called replica-symmetric (RS) cavity equation: 

1 

r? = g(r;,...,ry), 3(771, . . . , r,e_i)(a) = ■ (3) 

^Wl 1 ■ • ■ ; Vc~l ) T 
CTl ,...,CTe-l 1 = 1 

Here z{{rii}) is a normalization factor ensuring that 9{'n){'^) = 1- The cavity method then proposes an expression 
for the free energy density / in terms of the solution rj of this self-consistent equation. 

The equation Q and the associated RS prediction for the free energy density are correct when the spatial decay 
of correlations between variables is fast enough, which is equivalent to the existence of a single pure state in the 
Gibbs-Boltzmann measure (or of a finite number of them, related by simple symmetries). At low enough temperature 
and in presence of many constraints between variables (i.e. when c is large enough) this assumption breaks down. The 
Gibbs-Boltzmann measure is then split in an exponentially large number of pure states (or clusters), the correlation 
decay assumption underlying the RS computation is then valid only inside one pure state, not in the complete 
Gibbs-Boltzmann measure. The 1-step of replica-symmetry breaking (IRSB) cavity method deals with this situation 
by making a further ansatz on the structure of the pure states. One assumes that the typical number of pure 
states of internal free energy density / is exponentially large, at the leading order e^^^^\ with a rate function S(/) 
called complexity or configurational entropy. This function is computed via its Legendre transform [s^ (j){m,) = 
minj{/ — TI](/)/m}, the parameter m conjugated to / being called the Parisi replica-symmetry breaking parameter. 
The thermodynamic potential 0(m) is obtained from the solution of the IRSB cavity equation that generalizes ([3]), 
the order parameter solution of this self-consistent equation becoming P(77), a distribution of the RS probabilities rj 
with respect to the choice of the pure states weighted proportionally to e~^'^^^ : 

= | IX{'^P{vi)5{'n-g{m,---,Vc-i)) z{in,...,Tic-ir , (4) 

where Z is a normalization constant and the functions g and z are defined in Eq. One can check that the 

replica-symmetric equation ([3]) is recovered if P is a delta-peaked distribution. 

The structure of the Gibbs measure and its decomposition into pure states can be probed by the inter- and intra- 
state overlaps. The latter are defined, respectively, by: 



J2 1 AP{r,)AP{r,')^{a)Ti'{a) , =11 j dPiv)vi<y)vi^) , (5) 



where P is defined by the right hand side of Q, but with a product over c terms instead of c — 1. The energy function 
([1]) is symmetric under the exchange of two colors; hence for all a, J dP{'ri)'ri{a) = — !— , and qa is always equal to 1/q, 



col- 

Moreover, on the RS solution, this symmetry enforces ry(cr) = l/<Zcoi for all cr. The appearance of a non-trivial IRSB 
solution to the equation (j?]) can then be detected by the fact that the intra-state overlap qi takes a value strictly 
larger than 1/gcoi- 

Depending on the values of the temperature and of the degree c the pure state decomposition of the Gibbs-Boltzmann 
measure exhibits qualitatively different properties, that can be read on the type of solutions of the RS (O and IRSB 
cavity equations: 

• if the IRSB equation with m — 1 admits only a RS solution, then almost all configurations are part of a single 
pure state, the RS hypothesis is correct, and its prediction for the free energy density is valid. We shall call this 
situation a RS phase, or paramagnetic phase (P) here. 
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FIG. 1: Classical phase diagram of the coloring problem for ^coi = 4 in the connectivity-temperature plane, see for 
the actual numerical values. The solid orange line separates the paramagnetic phase (P) from the dynamical paramagnot (dP) 
while the dashed blue lines marks the boundary of the genuine spin glass (SG) phase. The model is defined only for integer 
connectivites (indicated by squares and circles); lines serve as a guide to the eye. The connectivites Cd = 9 and ck, Cs = 10 are 
defined in the text. Pink arrows indicates the connectivities that will be studied in the quantum case. 

• if the IRSB equation with m = 1 admits a non-trivial solution, and if the min in the definition of 0(m = 1) is 
reached at /int with I](/int) > 0, then the system is in a dynamic IRSB (dlRSB) phase, or dynamic paramagnetic 
phase (dP) here. Almost all configurations belong to pure states of internal free energy density /int, and there 
are e^^^-'''"'^ such pure states. It turns out that the total free energy density /int — TI](/int) coincides with the 
RS prediction, yet the splitting of the configurations into pure states has drastic consequences on the dynamics 
of such a model, that becomes non-ergodic (on sub-exponential timcscales), hence the name of dynamic IRSB 
phase. 

• if the IRSB equation with m = 1 admits a non-trivial solution, but predicts a negative complexity S(/int(?7i = 
1)) < 0, then the Parisi breaking parameter m has to be set to the value uis G [0, 1] such that S(/int(™s)) — 0. 
In such a IRSB phase almost all configurations belong to a sub-exponential number of pure states of internal 
free energy density /int(?7is); as the complexity vanishes this value is also the IRSB prediction for the total free 
energy density of the model. Thise phase corresponds to a genuine spin glass (SG) phase. 

In general these three phases are encountered in this order upon lowering the temperature (for large enough degrees) 
or increasing the degree (for low enough temperature). The transition between the RS and dlRSB phase occurs at the 
dynamic temperature Td(c) (or at Cd(T)), and is not accompanied by any singularity in the free energy. The transition 
between the dlRSB phase and the IRSB phase, at which the Gibbs-Boltzmann measure condenses on the few largest 
existing clusters, occurs at Tk^c). This transition corresponds indeed to the Kauzmann transition in the Random 
First Order Theory (RFOT) of structural glasses 53, 54] , which is a second order transition from the thermodynamic 
point of view, even if the order parameter jumps discontinuously ("first order") when the transition is crossed. In this 
context the dynamic transition corresponds to the divergence of the relaxation time of the Mode Coupling Theory 
(MCT) of supercooled liquids [5a]. T his p attern of transitions has indeed been found for the coloring of random 
regular graphs with gcoi > 4 colors [H, [11], and is illustrated on Fig. [TJ The case gcoi = 3 is singular, there is then 
no intermediate dlRSB phase, the transition between the RS and the IRSB phase occurs continuously, via a local 
instability of the RS solution towards the IRSB one at a temperature Ti(c). In the context of optimization problems 
the zero temperature limit plays a particularly important role; the transition lines defined above end in this limit 
at the critical connectivities and ck- A further transition can be defined at zero temperature: the satisfiability 
threshold Cg is such that the ground state energy vanishes for c < Cg (in terms of the coloring problem the typical 
answer to the decision question is yes, there are proper colorings of such graphs) and is strictly positive for c > Cg. 

In the following we shall study the quantum version of this model for QcoI = 4 and c = 9 or c = 13; these two cases 
correspond to the situations illustrated by arrows on Fig. [T] Before turning to this quantum framework let us make 
a few remarks. The coloring problem can be studied for other ensembles of random graphs, notably for Erdos-Renyi 
random graphs in which the degree of vertices are Poisson random variables of average c. The transitions described 
above are also encountered in this case [i^S^, with the difference that c is now a continuous parameter. We restricted 
ourselves to the regular case for technical reasons: to deal with the fluctuating degrees in the Erdos-Renyi case one 
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has to introduce a more complicated order parameter, i.e. in the IRSB cavity treatment a distribution over the 
distributions P to deal with the fluctuations of the probabilities rj with respect both to the choice of the pure states 
and of the local connectivity. Although this complication is affordable in the classical case, it would become impossible 
to treat in the quantum setting. We also neglected in this brief presentation the phenomenon of full replica-symmetry 
breaking 44]: the IRSB description of the Gibbs-Boltzmann measure is only the first one in a hierarchical construction 
in which the pure states are themselves grouped in clusters of pure states, and so on and so forth, that become relevant 
at a so-called Gardner transition , computed in [is", "4^ for the classical antiferromagnetic Potts model. This whole 
hierarchy can be dealt with in fully-connected models, most notably the Sherrington-Kirkpatrick one, yet it cannot 
be handled with present techniques in finite connectivity models (defined on sparse random graphs) , even for classical 
variables, not to speak about quantum models. The full RSB structure allows to study some important physical 
properties of these models, most notably the marginal stability of glass states; however, the IRSB approximation is 
expected to give quantitatively good predictions for thermodynamic properties even when it is unstable towards a full 
RSB one, and therefore it will be sufficient for the purposes of this study. 



B. Definition of the quantum version of the coloring problem 

In order to define the quantum version of the problem we introduce the g^j-dimensional Hilbert space spanned by 
G {1, • ■ • , 'J'coi}^}- We then define the "problem" Hamiltonian Hp corresponding to the cost function (P) as 
the operator diagonal in this (so-called computational) basis with diagonal elements equal to the classical energies: 

Hp = J2E{a)\a){a\ . (6) 

We introduce quantum fluctuations, i.e. off-diagonal matrix elements, by defining for each site i an operator Ti that 
flips the color ai to any other different color: 

W k') = '^"^^ ^"'^ = "'^ ^ V (7) 
I otherwise 

This is the simplest generalization from the Ising spin case (that corresponds to gcoi = 2) of the Pauli matrix af 
(taking as the computational basis the eigenvectors of the (if matrices) . The intensity of the quantum fluctuations is 
controlled by the "transverse field" P, the total Hamiltonian reading 

N 

H{r) ^Hp+THq, i/q = - ^f, . (8) 

1=1 

The partition function at inverse temperature /3 is Z(/3,P) = Tr[e~'^^'^''], thermodynamic averages being denoted 
(•) = Tr[« e-^"''^^]/Z{l3,T). In particular we will call rux = (Ti) G [O,qco\ — 1] the "transverse magnetization", by 
analogy with the Ising spin case. This also explains the name of "transverse field" for P. 



C. A sketch of the quantum cavity method 

A standard way to compute the partition function of such a quantum model consists in using the Lie-Suzuki- 
Trotter formula to disentangle the two non-commuting parts of the Hamiltonian i?(P), introducing copies of the 
original degrees of freedom ai along an "imaginary time" axis of length /3. In the limit where the number of such 
copies goes to infinity one thus obtains an exact representation of the quantum partition function as a path integral 
of a classical model, the price to be paid being the replacement of the spins ai € {!,... jQcoi} by more complicated 
degrees of freedom cr^ (we emphasize the use of a bold font here), which are piecewise constant periodic functions 
ai{t) from [0,^] to {1, . . . ,gcoi}- 

Apart from this replacement the classical model thus obtained from the quantum one has the same topology of 
interactions as the classical one ([T]); in particular if the latter can be treated with the classical cavity method (i.e. if it 
is defined on a locally tree-like random graph) , then its quantum version can be handled with the quantum extension 
of the cavity method where the spins ai are replaced by the imaginary time trajectories (Ti. This observation was 
first put to work in jJQ with a finite number of Suzuki- Trotter slices, the continuous imaginary time limit being taken 
in [42l |. These two works developed the quantum cavity method at the RS level, and the quantum IRSB framework 
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was then introduced in [2J, |57[. Detailed explanations and derivations, covering the case of the coloring in presence 
of a transverse field, can be found in the review [l7j . we shall thus content ourselves here with a few remarks, some 
technical details being presented in App.jBj in particular the parallelization of the numerical code we used to accelerate 
the resolution of the cavity equations. 

As explained above the quantum cavity method has exactly the same structure as the classical one, sketched in 
Sec. Ill Al The main modification is the fact that the basic object rj appearing in the RS Q and IRSB (U) cavity 
equations is now a probability distributions over the space of piecewise constant periodic functions from [0,/3] to 
{1, . . . , Qcoi}- This space is obviously much larger than {1, . . . , gcoi}, and in consequence a single 77, that could be 
represented by gcoi — 1 real numbers in the classical case, has now to be represented in an approximate way for the 
numerical resolution of Eqs. ([3]), (U). A convenient representation is provided by a finite sample of Atraj random 
trajectories, rj being approximated by the empirical distribution of this sample. 

Finally, the quantum overlaps are defined in analogy with the classical case, taking an extra average over the 
imaginary time: 



(9) 



qi = 

cr,cr 



As in the classical case, because of the symmetry between colors, qq = l/gcoi, and qi > l/^coi, with equality on the 
RS solution. 



D. Expected features for the quantum model 



1. Quantum phase transitions 



The ground state of the classical coloring problem is diagonal in the computational basis, while on the other side 
of the interpolation the ground state of the quantum Hamiltonian Hq is diagonal in the tensorial product of the 
eigenvectors of the fiipping operators T^. This defines two phases with very different physical properties. Slowly 
interpolating from high T down to F = 0, the spins have to rotate to go from the quantum paramagnetic phase to the 
classical spin glass one. If the classical Hamiltonian was simple or resembled the quantum one, for example if it acted 
as a product of single spin operators, this rotation would be easy and nothing special would happen along the way; 
but in our case, because of the complicated nature of the interaction © and of its associated classical phase, one can 
expect collective effects to become strongly relevant and a quantum phase transition Iss!] to occur along the way. 

It is well established that the gap vanishes at least polynomially fast with at a quantum second order critical 
point [2^, [13 (it vanishes exponentially fast in some cases in presence of disorder [s^), while it vanishes exponentially 
in at a first order phase transition |23 - l24| . As the scaling of the minimal gap is, according to the adiabatic theorem 
, critical for the success or failure of the QAA, we have to determine in our case the order of this quantum phase 
transition. This will be the subject of Sec. IIIII 



2. Level crossings and the role of entropy 



Apart from a quantum phase transition at relatively large values of the transverse field, we also expect the addition 
of the transverse field to have very important effects also for much smaller F. A first example of such a phenomenon 
was put forward in [30l l3lj |: we sketch it very briefiy here. Consider an instance of an optimization problem which 
has an isolated solution, and another local minimum of its cost function with only one violated clause, and such that 
these two configurations are far apart in Hamming distance. Computing the continuations of these energies using 
standard perturbation theories, one can find that, if the model is well chosen, they cross for some value of F. This 
value shrinks to zero as A^ — > 00, hence the name of perturbative crossings in the literature. This mechanism is for 
the moment only understood in perturbation theory, and therefore it holds whenever perturbation theory holds, that 
is, if at small enough F the full eigenstates of the quantum problem remain close enough to the classical eigenstates. 
The other important ingredient is a very particular construction of the instances of the problem, that admit only one 
solution. But typical instances of generic random optimization problems, even close to the satisfiability threshold, 
have an exponentially large number of solutions, and so we expect that non-degenerate perturbation theory should 
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not hold [601 and that the spectrum should be much more complex. Therefore we would like to understand what 
happens generically in problems that have multiple and not necessarily isolated solutions. 

The quantum random subcubes model [TtI . l32j is built as a quantum extension of the classical random subcubes 
model [61] (note that this model makes use of Ising spins, but these could be changed into Potts spins with any 
number of states gcoi, with only irrelevant changes in numerical prefactors). In this model, the classical spin glass 
phase is represented by a set of random subcubes of the Hilbert space (representing clusters). These subcubes are 
supposed to be disjoint, of different sizes, and to be associated with random (classical) energies. Adding quantum 
fluctuations to this model gives rise, under reasonable hypotheses on the distribution of the sizes of the subcubes, 
to a series of level crossings induced by a combined energetic-entropic effect. More precisely, if one introduces the 
internal intensive entropy s of a cluster (such that a cluster contains e^" configurations), the complexity E(e, s) such 
that there exist e-^^^^'") clusters of intensive energy e with such an entropy, and Sinax(e) = sup{s, E(e, s) > 0}, then 
it can be shown that the ground state energy of the model is given by: 

ecs (r) = min [e - r(ln 2) w(e)] (10) 

e 

Henceforth, as soon as F > l/((ln2)s,'„a^x(cGs(r = 0))), the minimum is obtained for a different value of e for each value 
of r. In this region the ground state changes abruptly from one cluster to another upon changing F by an infinitesimal 
amount, similarly to what is called temperature chaos in classical spin glasses [62, [63|- Because the clusters are at 
Hamming distances proportional to A^, we expect these crossings to be avoided at finite N by producing exponentially 
small gaps. 

The analysis of the quantum random subcubes model can also be done at finite temperature, which reveals the 
existence of a condensation transition at a temperature Tk^T) (similar to the one discussed in Sec. Ill A]l . The effect 
of the quantum fluctuations is particularly drastic when the classical model, at F = 0, remains un-condensed down to 
zero temperature (i.e. when rK(F = 0) = 0). Indeed in this case, for infinitesimal values of F > the model suddenly 
condenses on the largest clusters of classical groundstates, because the latter undergo the largest decrease in energy 
when the transverse field is turned on. One can study more precisely the limit F,T — !■ 0, in which the partition 
function is equivalent to 



ds exp 



^Mvr (V m ln(2cosh(/3F)) ^ 
N [ S(eGs(r = 0), s) + s — I 



(11) 



The slope of the condensation line 2k(F) in the small F limit can then be deduced from a saddle-point evaluation of 
this expression: it corresponds to the critical value of /3F such that the maximum of the argument of the exponential 
is reached at the upper limit of integration Smax(eGs(r ~ 0)), i.e. 

Tk{T) ^^^F/Argcosh[2-^=^(<=«s(r=")^'^— (12) 

As discussed in Sec. Ill A[ we expect the coloring problem for q^oi = 4 and c > 9 to possess a clustered classical 
phase and thus to exhibit a phenomenology very close to that of the quantum random subcubes model; in particular, 
crossings within the spin glass phase (see Sec. IIV[) and a linear condensation temperature Tk(F) at small field for 
c — 9 (see Sec. IIII B|) . The fact that we shall observe these phenomena in the coloring problem confirms the relevance 
of this scenario, first proposed in [TtI. Is^. for generic random optimization problems. Let us finally briefly explain 
why the same properties did not appear when studying the random regular Xorsat problem |24i] , even when it has a 
degenerate ground state (for example 4-Xorsat on a regular random graph of connectivity 3): because such a formula 
is locked [6J|, each solution is isolated and no entropic effects appear when adding quantum fluctuations. Hence the 
quantum coloring case is the first realistic optimization problem displaying the phenomena outlined above. 



III. THE PHASE DIAGRAMS OF THE COLORING MODEL IN A TRANSVERSE FIELD 

In this section we present our results concerning the phase diagrams of the quantum model as a function of the 
temperature T and of the transverse field F, for two different connectivities c, emphasizing the modifications of the 
classical transitions (at F = 0) once the quantum fluctuations are turned on. Most technical and numerical details 
of the analysis are deferred to the Appendices. A further investigation of the structure of the spin glass phase is 
presented in Sec. |IVl 
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FIG. 2: Phase diagram of the coloring problem for q^oi = 4 and c = 13. The solid black line is the continuous transition line 
between the IRSB spin glass phase and the RS (quantum) paramagnetic phase. For numerical accuracy reasons we were not 
able to determine the continuation of the dynamical Td and Kauzmann Tk transitions when F > 0. 

A. The Qcoi = 4, c = 13 case 

Let us start with the case c — 13. We recall that in the classical case (F = 0) the model exhibits a clustering 
(dynamical) transition at rd(r = 0) = 0.455, a condensation transition at 7k (r = 0) = 0.450, and a RS instability 
transition at Ti(F = 0) = 0.441 [H] (we will not discuss the Gardner transition towards full RSB at Tq = 0.185 jSj). 
From a thermodynamic point of view the only relevant transition in this classical limit is Tk, the free energy has a 
discontinuity in its second derivatives at this temperature, while it is non-singular at and Ti. The latter marks the 
limit of local stability of the RS solution, but it is here irrelevant because of the discontinuous transition towards the 
IRSB phase that occurs at higher temperature. 

We report on Fig. [2] the phase diagram that we obtain with the quantum IRSB cavity method. We do find, for 
F large enough, only one transition line, 71 (r), that corresponds to a local instability of the RS solution towards a 
non-trivial solution of the IRSB equation. This transition is continuous, in the sense that the overlap qi defined 
in Eq. ([9]) grows continuously from its replica-symmetric value 1/gcoi when one enters the spin glass IRSB phase. 
In order to locate it as precisely as possible we resorted to a finite population size scaling analysis, as explained 
in more details in Appendix IC II On this line the free energy of the model is singular, and this corresponds to a 
thermodynamically third order transition. This pattern of a second order RFOT like condensation transition (at Tk 
for F = 0) that becomes a thermodynamically third order transition upon changing one parameter was actually met in 
earlier studies of mean-field disordered models, namely in the fully-connected spherical p-spin model in presence of a 
magnetic field [H, [6^. For the convenience of the reader we gathered in Appendix [X\ the main aspects of the analysis 
of this (classical) model; the third order character of the thermodynamic transition derives, in both cases, from the 
linear growth of the overlap with respect to the distance to the continuous transition. A further evidence for the order 
of the transition is presented on Fig. [31 on the left panel one observes a rather good collapse of the complexity curves 
S(to), for various values of the transverse field F at a fixed temperature T, rescaled by (TiiT) — F)^. These curves 
vanish at the static Parisi parameter ms(T, T), that admits a finite limit (smaller than 1) when F — > Ti{T). As shown 
on the right panel of Fig. [3] this value is proportional to T when T is reduced for a fixed value of F. 

In the classical limit F — 0, the line Ti(F) falls below the classical transitions Td and Tk- This means that these 
two classical transitions should extend into two lines at finite but small F, and at some finite F the transition must 
change nature from discontinuous to continuous. However, in this model all this happens in a range of T and F too 
small to be resolved within the accuracy of the numerical resolution of the quantum IRSB equation. 

We should emphasize here that the incorporation of the replica-symmetry breaking effects was crucial to unveil the 
phase diagram of this model. As a matter of fact the RS computation predicts, incorrectly, a first order phase transition 
as a function of the transverse field at low enough temperature. We present in Appendix IC 31 the detailed numerical 
results and arguments we used to rule out this spurious prediction of the RS computation. A further confirmation on 
the continuous nature of the RSB transition will be provided by the Monte-Carlo simulations presented in Sec. IIVI 
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FIG. 3: Results for the coloring problem for gcoi = 4 and c = 13. Left panel: scaling form for the complexity upon approaching 
the third order transition. The temperature is fixed to T = 0.06, while F is varied. Upon approaching the critical field 
Fi(r = 0.06) = 0.595, the curves collapse. Note that the static value of m, ms(r, F), at which the complexity E(r, F,m) 
vanishes, is not singular in F close to the transition. 

Right panel: static value of the Parisi parameter divided by T, for various values of T. The dashed vertical line is the zero 
temperature limit for the critical field, Fi(r = 0) = 0.598. In the low temperature limit, curves collapse. The black (resp. red) 
dashed lines are fit to ms(T, F)/r for T = 0.06 (resp. T = 0.09) near the critical field. 

B. The qcoi = 4, c = 9 case 

We now turn to the case c — 9. Before discussing the quantum case, we recall ^] that for gcoi ~ 4 and c = 9, the 
model is classically satisfiable, meaning that the ground states' energy is zero and that graphs are typically colorable. 
These ground states are exponentially numerous and are arranged in an exponentially large number of clusters; this 
corresponds to the region Cd < c < ck, meaning that TK(r = 0) = 0, while Td(r = 0) = 0.153 is finite. The model is 
thus described, for F = and T < Td, by the IRSB equation with Parisi parameter m = 1. 

1. Solutions of the cavity equations and spinodal lines 

As discussed in details in [13, 113 1 the central quantity that is computed by the cavity method is the free energy 
function (j){ni), which is the Legendre transform of the complexity function, 0(m) = miny{/ — TT,(f)/m}. The 
thermodynamic free energy of the system is obtained by maximizing (j){m) in the interval ni G [0, 1], therefore in the 
following we restrict ourselves to these values of m. 

When solving the cavity equation, we found that in addition to the trivial RS solution, in some regions of (T, F) there 
exist two different non-trivial IRSB solutions. The first one develops continuously from the RS solution through a 
linear instability: hence the overlap qi grows continuously from its minimal value 1/gcoi, and we refer to this solution 
as the "low overlap" solution IRSBiq. The second solution appears discontinuously, as in the classical coloring 
problem: the overlap qi is typically larger, hence we refer to this solution as the "high overlap" solution IRSBhq. 
The coexistence of two different IRSB solutions is quite unusual, but had been observed before: a good pedagogical 
example is discussed in details in Appendix \K\ 

We find that the IRSBiq solution exists in an interval m S [0,771^^], while the IRSBhq solution exists in an interval 
m G [jT^hq' 1]' illustrated in Fig.|4l Both rn^ and to^^ depend on T and F, to lighten the notations we keep implicit 
these dependencies. To delimit the region of existence of these solutions in the (T, F) phase diagram, we define the 
following transition lines: 

• The clustering (or dynamic) transition line Td(F) is defined by the condition m^^ = 1. It is the continuation of 
the classical clustering transition Td(F = 0) to the quantum regime. Above this line, the IRSBhq solution only 
exists outside the interval m € [0, 1] and is therefore irrelevant for the thermodynamics of the model. 

• The transition line Ti(F) is defined as the point where the RS solution becomes linearly unstable. Note that 
the instability of the RS solution is independent of m and leads to the IRSBiq solution, hence we cannot have 
coexistence of the RS and IRSBiq solutions. Above the line Ti(F), the IRSBiq solution does not exist as it 
coincides with the RS one. Below this line, the RS solution does not exist as it becomes the IRSBiq one, at 
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FIG. 4; (Left) Spinodal lines of the coloring problem for q^oi = 4 and c = 9. The orange line with squares is the dynamic 
transition line Td(r), the black line with ticks is the continuous transition line ri(r), the red line with crosses is the line T*(r). 
Lines are fits to the data which are shown as guides to the eye. (Right) Sketches of the function ^(m) corresponding to the RS 
and the two IRSB solutions in the three regions defined by the intersections of the lines in the left panel. For instance, this 
would correspond to fix T = 0.1 and reduce F from top to bottom. 



least in the interval m e [0, ml^]. The temperature Ti(r) is non-monotonous, hence its inverse function has two 
branches that we denote by rj'"(r) and rr(T). 

• A third line Ti,{r) is defined by the condition mj^ — m^^. When this happens, the two distinct IRSB solutions 
merge into a unique IRSB solution. Note that (f>{m) is a convex function of m. When mj'q = m^^ ~ rrii, the two 
solutions merge continuously and with continuous first derivative, (/)iq(mj,) = 0hq(w*) and 4>[^{mi,) — (f>'^^{mi,). 
In this sense the line Ti,{T) corresponds to a "critical point": below this line, there is a first order transition 
in m between the two non-trivial solutions; on the line, the first order transition disappears, a transition still 
exists and is of second order; above the line, there is a unique analytic IRSB solution that we still call IRSBiq 
because it connects continuously to the RS one. 

As illustrated in Fig.|4j these three lines divide the (T, T) plane in four different regions. At high enough temperature, 
T > max(Ti, Td), only the RS solution exists. In the region below and above Ti, the RS solution remains stable but 
it coexists with the IRSBhq one (right, lower panel). In the region below min(rd, Tj, T,,), the two non-trivial IRSB 
solution coexist (right, middle panel). In the region below Ti and above max(Td,T-t), only the IRSBiq exists, either 
because the IRSBhq has moved outside the interval m e [0, 1] at Td, or because the two solutions have merged at T^, 
(right, top panel). Note that in the numerical computation each solution, for the same point of the phase diagram, 
can be reached by following different routes: starting from high T and decreasing it to reach the IRSBiq solution, or 
starting from T = and increasing it to reach the IRSBhq solution. 



2. Thermodynamic phase diagram 

The thermodynamic free energy of the problem corresponds to the maximum of 4>{m), which is reached in m = 
ms(T,T). We have seen that in some regions of the phase diagram several solutions for (j)(m) can coexist, leading 
to different branches of this function. This leads to several thermodynamic phase transitions that we now describe. 
For clarity, in Fig. [S] we report in the left panel the thermodynamic lines only, while in the right panel we report the 
complete phase diagram including the spinodal lines that were already discussed in Fig. 01 

Let us recall the following properties of the IRSB solution. First of all, 0(m = 1) is always equal to the replica- 
symmetric free energy. Second, the complexity T,{m) = I3m?(j)'(m) is related to the first derivative of (f>{m). We find 
the following phase transitions in the model: 

• At high temperature the system is in a RS paramagnetic (P) phase. Upon lowering the temperature at low 
enough T, the dynamical transition Td(r) is met and below it the IRSBhq appears. As in the classical limit 
r = 0, for r small enough the new solution is such that S(m = 1) c>c 0'(m — 1) > 0. By convexity, for any 
m < 1 the free energy of the IRSBhq solution is smaller than the one of the RS solution and the system remains 
in the paramagnetic phase. However, as in the classical case the paramagnetic phase is not ergodic because it 
is a superposition of many metastable states, and we refer to it as "dynamical paramagnet" or dP. 
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FIG. 5: Phase diagram of the coloring problem for gcoi ~ 4 and c = 9. In the left panel only the thermodynamic lines and 
the thermodynamic phases are reported. The right panel contains all the lines reported in the left panel as well as in Fig. [J] 
when two solutions of the cavity equations coexist, the stable one is underlined. The orange line with squares is the dynamic 
transition line rd(r), the black line with ticks is the continuous transition line T'i(r), the red line with crosses is the line 7+(r), 
the blue line with circles is the condensation transition Tk(r). Lines are fits to the data which are shown as guides to the eye. 
In addition, we report a purely conjectural shape of the first order transition T'fo(r) as a purple full line. At large F, ri(F) ends 
at F+(r = 0) ~ 0.47. 

• Upon lowering further the temperature, the complexity S(m = 1) decreases and reaches zero at the condensation 
temperature TK(r). At this point the slope of (jibqim) in 1 becomes negative, in such a way that 0hq('7i) has 
a maximum in the interval m S [0, 1], and its value at the maximum is necessarily larger than the RS solution 
(Fig. m right, lower panel). Hence at Tk a thermodynamic transition happens between the RS and IRSBhq 
phase. This transition is a standard RFOT transition: it is of second order thermodynamically, because the 
maximum of (t>hq(jn) moves smoothly enough away from to = 1, but it is of first order from the point of view 
of the order parameter because the system is jumping discontinuously from the RS to the IRSBhq phase. We 
call the spin glass phase below Tk the SGhq phase. Note that we find numerically that the line T'K(r) grows 
linearly with F at small F: this is a consequence of the clustered structure of the classical spin glass phase, as 
will be discussed in Sec. IIV Al 

• If the temperature is lowered below Ti(F) at high enough F, the RS phase becomes linearly unstable before the 
Une Tk is met. In this case, the possible presence of the IRSBhq solution is irrelevant as we already discussed. 
On the other hand, the solution IRSBiq appears at Ti(F) and we always find that it appears with a negative 
complexity at to = 1, I](to = 1) < 0, which implies that the slope of 4)iq{'m) in m = 1 is negative. The 
maximum of (f>{m) is thus assumed on the IRSBiq in this case (Fig. right, upper panel). The transition at 
Ti(F) is a thermodynamical third order transition (see Appendix . and the overlap grows continuously from 
its replica-symmetric value. We call the phase below Ti the SGiq phase. 

• Finally we should discuss the behavior in the spin glass region, for T < min(TK,Ti). The line Tk(F) can be 
continued in this region, and above it 0hq(TO = 1) > 0, hence the IRSBhq is surely metastable with respect to the 
IRSBiq solution. However, below this line, 0hq(TO = 1) < 0, hence both the IRSBiq and IRSBhq solutions have 
a maximum for m e [0, 1]. This can lead to a thermodynamical first order transition between these two solutions 
if the values of (/)(m) at the two maxima cross. This must happen at some temperature 7fo(F) < Tk(F), because 
we know that the IRSBiq solution is stable at high temperature while the IRSBhq is stable at low temperatures. 
Unfortunately, the free energy differences are so small in this region that our numerical accuracy does not allow 
us to determine the first order transition line. The line drawn in Fig. [S] is therefore schematic. We can only say 
that the variation of the complexity at m = 1 is found to be much higher in the IRSBhq solution than in the 
IRSBiq one. This is actually related to the fact that Tk is a second order transition while Ti is a third order 
transition. It implies in particular that Tfo must be tangent to Tk when they separate, and it also suggests that 
Tfo(F) must be always very close to Tk(F). Note also that the line Tfo(F) must necessarily end on the line Tj,(F) 
where the distinction between the two IRSB phases disappears. 

As a final remark, we note that within our numerical accuracy it seems that the three lines Td(F), Tk(F) and T*(F) 
cross at a single point. This finding is perfectly compatible with all of our data. To better understand this point let 
us call (To,Fo) the point where the lines Td(F) and T*(F) cross. If Tk(F) does not cross the lines in the same point. 



12 



there are two possibilities: (i) Tk crosses T^, at some T < Tq. However, this can be excluded because at the two 
IRSB solutions must merge and convexity implies that the slope of 0hq in m = 1 must be negative, (ii) Tk crosses 
Td at r < Fq. This scenario seems to us logically possible. However one must stress that in most cases the lines Tk 
and Td do not cross, but merge into some kind of critical point. Therefore the scenario of a triple crossing at (To, To) 
seems likely. Unfortunately, we could not devise a more solid argument. 

In summary, the thermodynamic phase diagram contains a paramagnetic (P) phase, a dynamical paramagnet (dP), 
a "high overlap" spin glass (SGhq) and a "low overlap" spin glass (SGiq). The transition at Td from the P to the dP 
phase is a standard dynamical (clustering) transition, and has no thermodynamic consequences. The transition Tk 
from dP to SGhq is also a standard RFOT transition, i.e. a second order transition. The transition Tj from P to SGiq 
is a third order thermodynamic transition. Finally, a first order thermodynamical transition between SGiq and SGhq 
must exist, even if our numerical accuracy is not sufficient to determine it precisely. 

In particular we find that in the limit of zero temperature, the system becomes a spin glass under the action of an 
infinitesimal transverse field. At very low temperatures the two spin glass phases transform smoothly into each other, 
and the only phase transition is the third order transition at T'^{T — > 0) ~ 0.47. 

The artifact of the spurious first order transition predicted by the RS computation, discussed for c = 13 at the end 
of Sec. IIII Al and in Appendix IC 3[ occurs also in the case c = 9, even if at even lower temperatures with respect to 
c= 13. 



IV. STRUCTURE OF THE SPIN GLASS PHASE 
A. Clusters crossings in the spin glass phase 

As sketched in Sec. lIID~2l apart from the various transition lines that appear for this model when the quantum field 
is turned on, we also expect the spin glass phase to exhibit an interesting quantum behaviour due to its rich classical 
structure. The intuition based on the study of the quantum random subcubes model [13, [111 recalled in Sec. Ill D 21 
is that a continuum of (avoided) level crossings occurs in the spin-glass phase, because the classical model admits a 
pure state decomposition with clusters of various energies and entropies. As clusters with a larger classical energy 
have also a larger entropy, their total (quantum) energy is reduced faster when F increases, hence the crossings which 
lead to exponentially small gaps because of the extensive Hamming distance between clusters. The linear behavior of 
the hue Tk(F) reported in Sec HlTBl was a first confirmation of this mechanism, we report here a further evidence in 
its favour (the following discussion was already partially undertaken in [TtI). 

The coloring problem for gcoi = 4, c = 9 is in its clustered phase at zero temperature and zero field. Its clusters have 
internal entropies that are distributed according to a large deviation function (the complexity), Af{s) = cxp[iVI](s)]. 
Typical configurations are found in clusters that have zero energy and value of the entropy s* such that I]'(s*) = —1, 
and the corresponding complexity E* = S(s*) is strictly positive. In particular, for this model E* ~ 0.012 and 
s* ~ 0.08. However, many other clusters with larger and smaller entropies exist, as well as clusters with positive 
energies. When F > 0, each cluster A of degenerate classical states transforms continuously into a set of quantum 
states, the lowest of which (the "ground state of cluster A" \GS{A))) has an energy per spin 

e(A,F) = eei(A,F)-Fm,(A,F) , 
eci(A, F) = {GS{A)\Hp\GS{A))/N , ^^3^ 
m,(A,F) = {GS{A)\Y,nGSiA))/N . 

i 

The key observation is that, as in the quantum random subcubes model [3^ . we expect mx{A,r) to be finite when 
F — >■ 0, \imr^omx{A,r) — m'^{A), because there exist degenerate classical ground states at Hamming distance 1 
inside A, and we expect m'^{A) to be positively correlated with the classical entropy of the cluster. At the same time, 
eci(A,F) cx F^ at small F. 

Therefore, the energy of the ground state of a cluster is linear at small F, e(F) ^ — FmJ^(A). Largest clusters yield 
the greatest decrease in energy when quantum fluctuations are switched on, and they dominate at zero temperature 
as soon as F > 0. Because these are the states with maximal entropy, they correspond to E(smax) — 0. Hence as soon 
as F > 0, the zero temperature complexity abruptly drops to zero. In other words, we expect the system to condense 
into the largest clusters under an infinitesimal amount of quantum fiuctuations. This in particular implies that, as in 
the quantum random subcubes model, a non-zero Tk(F) should emerge, and that Tk(F) oc F for small F, the effect of 
the transverse field entering the computation of the free energy under the combination /3F. This feature is confirmed 
by the phase diagram presented on Fig. [S] 
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FIG. 6: Transverse magnetization obtained by path-integral Monte-Carlo simulations with = 1000 spins. All the simulations 
were started from the same planted state at T = 0, and then runned increasing F at fixed T. The black line is the extrapolation 
at T = 0, showing that rni{A) = limr->o limT->o rn^iA, T, F) > 0. 

Let us now probe the hypothesis mJJ(A) > 0, first numerically. Our path integral formalism does not allow to work 
directly at zero temperature, hence we had to run several Monte-Carlo simulations at low temperature and extrapolate 
the results to the limit T — > 0. In this case, because the system is in a dynamical IRSB phase at F = 0, it is possible 
to use "quiet planting" [67] to construct configurations equilibrated at T = 0. We started all the simulations in the 
same quietly planted configuration at F = 0, for the same instance of the coloring problem; in this way we assumed 
that the simulations all follow the evolution of the same cluster. Extrapolating the results to T = gives the ground 
state properties of the cluster as m^{A) = limr^o liinT->o rnx{A, T, F), as shown on Fig. [51 Note that at the very low 
temperatures we investigated, and given the size of the graph, the planted solution at T = is actually an equilibrium 
configuration of the model at F = for all temperatures of Fig. [6) 

It is also possible to get a theoretical understanding of the existence of this finite limit for limr^o limT->.o nixiA, T, F). 
First of all in the quantum random subcubes model, one can derive the exact expression ( A, r,F) = ^tanh(F/T). 

Hence m^{A) = and a perturbation theory in F to compute mx{A,T,r) will be divergent when T — 0, the 

coefficient of the term linear in F being proportional to mx{A)/T. For the coloring problem, no closed expression can 
be derived for mx{T, F), but it is still possible to compute it perturbatively in F within the path-integral representation 
of the partition function. The details of the computation are given in Appendix [Dl The important conclusion is that 
in the limit of small F and T (with the limit F — > taken before T — > 0) one gets the asymptotic 

F 1 ^ 

mx{A, T,r)^--J2 (9auth(a, i) - 1)^ , (14) 

i=l 

where a is an uniformly chosen groundstate of the cluster A, and gauth(2j*) denotes the number of colors the site i 
can take without creating a monochromatic edge if the rest of a is kept fixed (in other words the number of colors 
that do not appear in the neighbors of i in a;). The expansion with F — before T — ^ is thus the same as for 
the quantum random subcubes model. Altough we cannot prove that m'^ (i.e. with the order of the limits reversed) 
is again proportional to ((?auth — 1)a in this case, we expect that it will still be positively correlated with it, and 
with the internal entropy of the cluster. This could be checked numerically by repeating the above computation for 
many different clusters of different entropies. In any case the demonstration given above of a non-zero m'^(A) is an 
additional piece of evidence in favour of the mechanism of crossings induced by the competition between the energy 
and the entropy of the classical pure states. 

B. Quantum Monte Carlo annealing 

Let us conclude this section by pointing out the consequences of this clustered structure of the spin glass phase for 
a quantum annealing of the problem. One expects that, because of this complicated structure, a quantum annealing 
of the problem starting from high F down to F = will remain stuck in large clusters of higher energy when F gets 
small, and will not be able to find the smaller clusters that contain the ground states. However, it is not possible to 
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FIG. 7: Energy (left panel) and transverse magnetization (right panel) as a function of F at fixed temperature T — 0.06 for 
Qcoi = 4 and c — 9. Quantum cavity computations are reported as blue solid lines, path-integral Monte-Carlo simulations (for 
a single system of size A'^ — 10000) as dashed black lines with triangles, both for increasing and decreasing F. The inset of 
the left panel shows a zoom on the region of low F where hysteresis is observed in the Monte-Carlo simulation. The red arrow 
indicates the value of the classical energy (0.0019) that corresponds to a classical annealing starting from Td. 

simulate the real-time Schrodinger quantum dynamics of such a diluted model for reasonable values of N (the size 
of the Hilbert space is equal to 4^). Still, the effects of clustering can be seen on a classical alternative to quantum 
annealing, namely the annealing of a path integral Monte-Carlo computation (PIMC), whose results are shown on 
Fig. [71 These simulations were performed in two ways: in the first run, we prepared a typical graph tog ether with one 
of its typical configurations at a very small temperature T — 0.06 via the quiet planting technique [67|; we initialized 
the PIMC at F = in this configuration (which, given the size of the graph, is actually at zero energy), and then slowly 
increased F. Note that the results of this run perfectly match the results of the IRSB quantum cavity computation, 
confirming the validity of our analysis. In the second run, we initialized and equilibrated the PIMC at F = 2 (in 
the paramagnetic phase), and then decreased F down to F = 0. On the scale of the figure, no difference between 
the two PIMC runs is observed, however a closer look (inset of left panel in Fig. [7]) reveals that decreasing F one 
obtains a positive residual energy at F = 0. The latter is found to be larger than the residual energy after an infinitely 
slow thermal annealing, which is obtained by preparing a quietly planted configuration at Td and performing a slow 
classical annealing down to T = (gtI , [68| . Although a direct comparison is not possible because the PIMC annealing 
has not been extrapolated to the infinitely slow limit (we used steps of transverse field dT = 0.02 with 1000 sweeps 
per step) , this result suggests that an annealing of a PIMC simulation is not more efficient than a thermal annealing 
for this model. We believe that this is once again related to entropic level crossings inside the spin glass phase, as in 
the quantum random subcubes model. Note that the residual energy of the PIMC annealing is a priori not directly 
related to the one of a true quantum annealing. 

Similar results are obtained for connectivity c = 13, however in this case the quiet planting technique cannot be 
used at T = because Tk > 0. Therefore we used as a starting point for the increasing F run a classical configuration 
obtained through a slow thermal annealing from Td. These results are reported on Fig. |S1 

V. CONCLUSIONS 

In this paper we presented a study of the classical coloring problem, to which a quantum transverse field is added 
to represent the action of a quantum computer performing a quantum annealing with the aim of finding solutions to 
the classical problem. The study of this particular problem was motivated by the fact that among classical random 
optimization problems, it is the simplest one that shows an exponential degeneracy of solutions inside clusters (in 
the jargon of [6J] , it is a non-locked problem) . Therefore, this is the simplest non-trivial model where the predictions 
of [32], that were obtained on a toy model, could be tested. We believe that similar conclusions would be reached on 
similar problems, such as random if-SAT. 

Our main results are the characterization of the phase diagram in the (T, F) plane of parameters, and a further 
description of the spin-glass phase at low temperature and transverse field. Concerning the latter issue, we found 
evidences that the mechanism described in [32| is at work also in the case of the coloring, more precisely that (i) 
at low F, the line Tk(F) is linear in F for c = 9, with Tk(F = 0) = and Td(F = 0) > (see Fig. [5]); and that 
(ii) the transverse magnetization of the ground state of a cluster goes to a positive constant in the limit F — ?> (see 
Fig. [5]). Both these results are direct consequences of the exponential degeneracy of solutions in the clusters |3^, and 
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FIG. 8: Energy (left panel) and transverse magnetization (right panel) as a function of F at fixed temperature T = 0.06 
for gcoi = 4 and c — 13. Quantum IRSB cavity computations are reported as blue solid lines, path-integral Monte-Carlo 
simulations (for a single system of size A'^ = 10000) as dashed black lines with triangles, both for increasing and decreasing F; 
those obtained starting from F = were obtained after a slow thermal annealing of the system from rd(F = 0) (in this case, 
because Tk(F — 0) is finite, it is not possible to use quiet planting at low temperature), hence the small discrepancy in the 
energy for F = between the cavity computation and the Monte Carlo results. The (erroneous) replica-symmetric results (for 
the magnetization only) are shown in red with dash dotted lines and circles; the absence of hysteresis around the transition 
(F ~ 0.56) confirms that the first-order transition predicted by the RS computation is spurious. 



they show that the general scenario proposed in using a toy model is at work in realistic random optimization 
problems. 

A new feature of the present study, which could not be expected based on the analogy with the random subcubes 
model, is that the quantum phase transition between the spin glass and quantum paramagnetic phases at low tem- 
perature is of third order thermodynamically. Although we do not know the scaling of the gap at such a transition, 
the analogy with the results of ,25] leads us to believe that the gap will be polynomial right at the quantum phase 
transition. On the other hand, based on the results of [111, we expect that the spin glass phase is characterized by a 
continuum of level crossings with an everywhere exponentially small gap. It would be extremely interesting to study 
numerically the gap in given instances of this model; exact diagonalization is impossible due to the rapid growth of 
the Hilbert space with N (as 4^), however this might be doable through Quantum Monte Carlo following psi. l69l|. 
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Appendix A: A reminder on the fully-connected p-spin spherical model in a longitudinal field 

1. Definition and phase diagram 

During the study of the quantum version of the coloring problem we have encountered several phenomena (in 
particular a third order continuous transition towards a spin glass phase, and multiple RSB solutions) that also 
appear in a much simpler, classical disordered model, namely the fully-connected p-spin spherical model in a field [65l |. 
Analytical computations, and in particular expansions around the transition line, can be performed explicitly in this 
model and have constituted a very useful guideline for the analysis of the quantum coloring model. For these reasons we 
briefly collect in this Appendix the main results on this well-known model that are enlightening with this application 
in mind and refer the reader to [H, |6^, [70l - [73j | for more details on this model. 

Its Hamiltonian is 

N 

i?((Tl, . . . , (Tat) = — ^ >/n,...,ipCTii . . . (Tip — (Ti , (Al) 

l<il<---<ip<Af i=l 
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where the quenched couphngs Ji-^^...^i^ are independent Gaussian random variables of zero mean and variance p!/2iV"r' 



and the degrees of freedom ai are continuous real variables subject to the spherical constraint = The 

partition function is defined as the integration of the Gibbs-Boltzmann factor e~^^ over the A^-dimensional sphere 
of radius y/N, and the associated free energy density concentrates in the thermodynamic limit around its quenched 
average. The latter can be computed with the help of the replica method and reads 



(f){l3,h)= sup (t)irsh{l3,h;m,qi,qQ) 

0<m<l 
0<<3o<gi<l 



(A2) 



where the variational function is 



1 f 8"^ 

^irsh{8, h;m,qi,qo) = -—n + log(2^) + y [(1 - qf) + m{qf - q^)] + il3hf{l - qi + m,{qi - go)) 



- ^ + m 1 log [1 -q^] + — log [1 - gi + m{qi - go)] 

1 - gi + m(gi - go) m m 



(A3) 



In this expression m is the Parisi parameter discussed in the main text, while go (resp. gi) is the typical overlap 
between two configurations in different (resp. in the same) pure states. Note that the IRSB potential 0(to) discussed 
in the main text corresponds to the maximization of this function with respect to gi and go, with m fixed. 

The phase diagram of this model, shown on Fig. [HI is obtained by solving the maximization problem defined in 
(IA2IA3I) : depending on the values of I3,h the supremum of 0iisb is either reached in the subspace of parameters 
Qo — Qi (this corresponds to a replica symmetric situation), or on a non-trivial IRSB solution with go < gi. The 
phase transition that separates these two regimes (and that reveals itself as a singularity in 0(/3, h)) as a function 
of the temperature, changes qualitatively depending on the value of h. For h smaller than he the transition is of 
the RFOT type, as described in Sec. Ill AI there exists a dynamical transition temperature T^{h) and a condensation 
(Kauzmann) transition temperature T^ih) < T^{h), the free energy undergoing a thermodynamic phase transition at 
the latter. On this line T^ih) the overlap order parameter is discontinuous, yet the value of the static Parisi parameter 
nis, i.e. the one maximizing (jA3|) . is ms(/i, Tk(/i)) = 1, which makes this discontinuous transition second order from 
a thermodynamic point of view. On the contrary for h > (but of course not too large) there is a single transition 
temperature Ti(h), at which the order parameter of the RSB phase grows continuously with ms{h,Ti{h)) < 1, and 
which is thermodynamically of third order. The three lines Td,TK and Tj meet at h = he- In the following we give 
some details on the derivation of these properties of the phase diagram, and also on some features that are not directly 
relevant for its thermodynamic properties (a spurious first order RS transition and the coexistence of different IRSB 
solutions at TO 7^ TOs). 
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FIG. 9: Phase diagram of the spherical p-spin model for p = 3 (all values of p > 3 are qualitatively similar). The dashed orange 
line corresponds to the dynamic transition, the solid blue one to the Kauzmann (condensation) transition line, and the dotted 
black one to the continuous transition line. SG refers to the IRSB spin glass phase, dP to the dynamical paramagnet (IRSB 
phase with m = 1), and P to the paramagnetic (RS) phase. The three lines meet at the point {he, Tc). The orange dash-dotted 
line is the limit of coexistence of the two RS solutions, whose free energies cross at a spurious first order phase transition not 
shown on this figure. 
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2. The RS solution(s) and its stability limit 



Let us define the replica-symmetric variational free energy, function of a single overlap q, by substituting qo — qi — q 
in (jASP : one then sees that this expression is independent of m, and is equal to 

</.rs(/3,/i;g) = -^|l + log(27r) + ^(l-gP) + (/3/i)2(l-g) + ^ + log(l-g)| . (A4) 

The stationary points of this function of q are solutions of the following equation, 

tpqP~i + ^phr^^^. (A5) 

Depending on the values of the parameters /3, h there exists either one or three solutions to the RS equation (jASp . In 
the latter case the intermediate one corresponds to a minimum of and can be discarded, while the two extreme ones 
compete to maximize (jj^s and thus cross at a first order transition line. The boundary of the domain of existence of 
multiple RS solutions in the {h, T) plane is more easily expressed parametrically, as a function of q; on this boundary 
one has the relation (|A5p and in addition 

The orange dash-dotted line in Fig. [3] represents this spinodal limit of existence of multiple RS solutions. The first 
order transition between these two solutions occurs on a line, not shown on the figure, that starts at some T > Q,h — Q 
and joins the spinodals at their cusp. We shall see that this transition is irrelevant thermodynamically, as it occurs 
inside the IRSB phase. 

It is easy to check that if g is a stationay point of 4>rs, then the first derivatives of (push with respect to qoj <Zi and m 
vanish in (go, 91, w) = (g, q, to), for all values of m, i.e. that (g, g, m) is a stationary point of ipi-c&h- Let us discuss more 
precisely its nature, assuming g is a local maximum of (jj-^s- All second derivatives of (?!>irsb with respect to go,gi,m 
which involve at least one derivative with respect to to vanish, hence one can concentrate on the Hessian of (^irsb with 
respect to go and gi, and in particular on its determinant. A short computation leads to 



det IIess(m, g) — 



9i3o dqadqi 



dqodqi 



to(1 — to) 5^^ri 



2/3 

qa=q ^ 



2^- {l-qf 



(A7) 



As TO e [0, 1] and d'^(f)is/dq^ < on a maximum of the RS potential, the sign of the determinant of the Hessian can 
only change when the last factor of this equation vanishes, i.e. when 

(A8) 

The black line Ti{h) (or, alternatively, the two branches hf{T)) in Fig.lHlcorresponds to the parametric representation, 
as a function of g, of the conjoint solution of (lASp and of the RS equation (jASp . It starts from (T = 0, h^{T = 0) = 0) 



(corresponding to g = 0), reaches a maximum in ( = v'^ — ^p-i ,hc = J 2pp-2 ) (when g = ^^), before 



hitting the zero temperature line again at /i; (T = 0) = y — - (for g = 1). 

Outside the region enclosed by Ti{h), i.e. at high temperature or field, the determinant of the Hessian is positive, 
corresponding to two negative eigenvalues. A local maximum of 4>i-s then corresponds to a local maximum of (/)irsb in 
the larger RSB subspace of parameters, and such a RS solution is at least locally correct. On the contrary once Ti{h) 
is crossed one of the eigenvalues of the Hessian becomes positive, there is one direction that increases 4>iish starting 
from the RS local maximum, which is thus locally unstable (for all values of the Parisi parameter in [0, 1]). We will 
come back in more details in Appendix I A 41 on the behavior of the model around the instability line Ti(h). Let us also 
mention as a technical detail the issue of the stability of the RS solutions in their domain of coexistence. The latter 
is traversed by the small field branch of the instability line, [T] ; in the intersection of the domain of coexistence 
of two RS solutions with the interior of the instability domain T < Ti{h) both RS solutions are locally unstable. 
When two RS solutions coexist with T > Ti{h) the one continuously connected to the high temperature regime is 
locally stable, while the other one is locally unstable. In any case even the locally stable solution is irrelevant as the 
non-trivial IRSB solution will have a larger free energy, as discussed below. 
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3. The RFOT-like transitions at small h 



We have seen above one mechanism for the appearance of rephca-symmetry breaking, namely the transformation 
of a local maximum q of (/>rs from a local maximum of (jjirsh with qq = qi = q to a saddle. Another possibility is the 
appearance of a local maximum of (/)irsb with qg qi, i.e. far away from the replica-symmetric subspace (hence the 
"discontinuous" character of such a transition). At variance with the local instability, which is independent of the 
value of m, the occurence of the discontinuous transition is m-dependent. We shall start by considering the most 
important case m = 1, and come back on this m-dependence in Appendix I A 51 

At m = 1 the stationary conditions of (/)irsb, ^^f^ = ^g^"** = imply that go is solution of the RS equation (|A5p . 
while qi > qo verifies: 

/(<Zi) = /(9o) , with fix) = ^pxP-' - . (A9) 

Only for some range of parameters the equation f{x) = f{qo) admits a non-trivial solution qi > go- On its limit of 
existence one has the additional condition f'{qi) = 0, i.e.: 

^p(p-l)gr' = 7r^- (AlO) 

The line T^{h), corresponding to this limit of existence and drawn in dashed orange on Fig. 1^1 can be obtained 
parametrically: from (|A10|) one obtains /3 as a function of gi, then go is expressed in terms of gi as a solution of 
/(go) = /(gi), and finally the field is obtained from (jASp . go being solution of this RS equation. The relevant interval 
for this parametrization is gi € [^^, ^zf]- At the lower limit the line Td merges with the local instability T\ at he] 

in this limit the two extrema of / merge, gi — > go and (jAlOl) indeed reduces to (jA8|) . The upper limit gi = 

corresponds to /i = 0, go = 0, and one can find the explicit expression T^ih = 0) = \J ^2(p^iy^ - 

The non-trivial IRSB solution just found at to = 1 might be the global maximum of (/)irsb, or not. To test this one 
has to compute the derivative of 0irsb with respect to to, in this point. As discussed in Sec. Ill Al this is proportional to 
the complexity; when the latter is positive to = 1 is indeed a maximum of (/)irsb, otherwise there will be a maximum 
at a value TOs < 1. The Kauzmann, or condensation, transition separates these two regimes, and corresponds to the 
line drawn in solid blue on Fig. [HI 



4. The third order transition at large h 



We now come back in more details on the behavior of the free energy around the RS local instability. We want 
to justify that the perturbations in the overlaps go and gi with respect to their common RS value are linear in the 
distance from the instability line, and that the first discontinuities appear in the third derivatives of the free energy 
4>{P, h) (hence the transition is thermodynamically of third order). In particular, the complexity takes a scaling form 
in {5T -\- Sh)^ when one approaches the instability line (see Fig. [TU)) . These last two facts have been used in the 
analysis of the quantum coloring problem. 

As explained in Appendix lA 2[ when crossing the instability line Ti(h) the RS extremum goes from a maximum 
to a saddle (in the larger IRSB space of parameters). By definition the determinant of the Hessian matrix written 
in Eq. (|A7p vanishes in (h^Tiih)). More precisely, one can diagonalize the matrix in this point and find that it has 

'-(1-to) 

TO 



one eigenvalue with eigenvector 



and one strictly negative eigenvalue associated to the eigenvector 



TO 
1 — TO , 

Consider now some parameters (/3, h) inside the unstable region, denote grs(/3, h) the associated solution of the RS 
equation, and expand the IRSB potential around this point in powers of (5+, (5_ as 

5(j) = (?!)irsb(/3, /i;go = <7rs(/3, h)-{l~ m)S+ + mS-,qi = qrs{l3,h) + mS+ + (1 - m)S-,m) - cjirsiP , h; q^^iP , h)) . (All) 

As gis is a RS solution, the expansion begins with second order terms. In the following we shall only need the 
coefficients of the terms of order (5^ and 6'\_, denoted a/2 and b/3 respectively, with 



to) 



^p{p-l)q,,{f3,hr-^ 



1 



il-qUP.h)y 



b ^ -— to(I-to) 



(1 _ 2m)^p{p -l){p- 2)grs(/3, hf-' 



(l-g,,(/3,/i))3 



(A12) 
(A13) 
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FIG. 10: Scaling form of the complexity h, m)/5h^ of the spherical 3-spin model for P = 1/0.3 and h = hf{l3) - 5h. The 
solid black line was obtained from the analytical expansion given in Eq. (|A14[I . 



Consider now that (/3, h) is at a small distance, call it e, of the instability line; a short reasoning reveals that the 
coefficients of the terms in 5^ and 5+5- are of order e, while all others are of order 1 in this limit (in fact, note that the 
perturbation 5+ is in the direction of the vanishing eigenvalue of the Hessian) . One then sees that the maximization 
of 5(1) over 5+ and 5- wiU lead to 5+ = 0(e), 6- — O(e^), and that at the lowest order 5(j) is simply given by the 
maximization of a5\/2 + b5\/'i. We thus obtain in this limit 



12/^ („,+ (! _2rn)(H^fi^)' 



^p(p-l)fc(/3,/i)P-2 ^ 



2^^" — (l-fe(/3,/,))2 



(A14) 



where we made some simplifications using Eq. (IA8|) verified by q, the value of the overlap at the point of the instability 
line approached in this limit. The term in square brackets is positive in the RS unstable region, and vanishes linearly 
on the instability line. As (j)^s{f3,h) is regular across this line this demonstrates the third order character of the 
transition. The dependence on m of the IRSB potential can then be easily studied, in particular its maximum is 
reached in fhs{q) — ^^^—^2q~~^- This behaviour is checked numerically by the scaling plots of the complexity, see 
Fig. 1101 We emphasize that as 5+ — 0(e) on the maximum of 0, the overlap order parameter qi — qo grows linearly 
away from the instability. 

It could seem at this point that the above study is valid along the two branches h^{T); this is however not the 
case, as further considerations reveal. The maximization over the IRSB overlaps must enforce the condition qi > qo, 
which translates here into 5+ > 0. For the local maximum of a5+/2 + b5+/3 to happen on this side of the origin one 
must have 6 < (we always have a > in the unstable regime). Assuming m e [0, 1], and after some simplifications, 
one sees this condition to be equivalent to fhs{q) + m(l — 2fhs{q)) > 0. Two cases can now be distinguished: 

• The high- field branch K^{T) corresponds, in the g-parametrized representation, to q S 1 • Then ms{q) e 

[0, 1], thus the condition qi > qg at the local maximum is respected for all m G [0, 1], and the maximum of 0(m) 
is reached in an acceptable value fhs{q) € [0, 1] of the Parisi parameter. 

• On the contrary for the low- field part of the instability line h[ {T), rhs{q) > 1. The condition qi > qo is only 
fulfilled for m G [0,m^q], where = 2ja"(g)-i • these values of m there exists a non-trivial solution of the 
IRSB equations that grows continuously from the RS one. But for m G [rttl^, 1] one sees that the maximum of 
the expansion would correspond to 5+ — >■ oo, in other words the maximum of the IRSB potential is far from the 
RS solution, and corresponds to the discontinuous solution. 



5. Coexistence of RSB solutions 



We have just seen that close to the low-field branch (T) < he of the instability line the locally unstable RS solution 
gives birth continuously to a non-trivial solution for m G [0, (T, h)]. On the other hand we have also shown that 
for m = 1 a discontinuous IRSB solution exists for T < Td(/i); actually the latter exists for m G [m^ {T, h),l], hence 
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S(m) 




FIG. 11: Study of the coexistence of two RSB solutions, for T — 0.55 and a field h — 0.166 in the narrow range [h^ (T) ~ 
0.1649, /i* (T) ~ 0.167] where the two RSB solutions have a common domain of existence m £ ["^hq'^^^ql- ^^o'' ^ ~ 0-55 
and h — 0.166 the continuous solution exists for m G [0,ml^], with mj^^ ~ 0.48, while the discontinuous one is defined for 
me K^.l], with m=P -0.34. 

Left panel: Complexity h, m) as a function of m. The discontinuous solution has the largest complexity in the domain of 
coexistence [m^^, m^^]. Note that the complexity of the continuous solution is very small, but non-zero. The relevant value of 
the parameter m for the thermodynamics, ms{T,h) ~ 0.77 where the complexity vanishes, is outside the domain of existence 
of the continuous solution. 

Right panel: The potential V{m,qi) defined in HA15P as a function of gi for several values of m. The inset is a zoom near 
gi — qo, and the ticks indicate the replica-symmetric solution. 

two IRSB solutions coexist for m e ["^hqC^' "^iq ^)] (whenever m^^{T, h) < (T, h)). This leads the physical 
observables to have two branches when plotted as a function of m (see left panel of Fig. [TT|) . However, once one also 
extremizes over m, no sign of the continuous solution to the IRSB equation remains, and the physical solution is 
always the discontinuous one. To justify this point it is convenient to consider the reduced potential: 

V{m,qi)^ sup (/)ii.sb(/3, gi,go) - sup (I)ii-sh{l3,h;m,q,q) , (A15) 
o<«o<9i <je[o,i] 

where the dependencies on T and h are kept understood, and the normalization has been chosen such that V vanishes 
on the replica-symmetric solution. A plot of 1^ as a function of qi for h slightly larger than h^{T) and several 
values of m is shown on the right panel of Fig. 1111 It can be seen that, as expected from the previous discussion, 
V has for m small enough (m < m^^{T,h)) a local maximum near the replica-symmetric point qi ~ qo (indicated 
by a tick on the figure), corresponding to the continuous solution. On the other hand, there exists for m sufficiently 
large (m > ni^^(T,h)) another local maximum of at a value of qi which is further away from go- This second 
maximum corresponds to the discontinuous solution. When m — )■ (T, /i), the continuous solution merges onto the 
discontinuous one. Increasing h, m^^(T, h) shrinks and the limit of coexistence of both solutions (at fixed m) is given by 
the critical field K{T) such that m^^{T,K{T)) = (T, /i^(r)). In the region h G [hr {T),hi,{T)] thermodynamical 
observables acquire two branches when computed at fixed m, as shown for the complexity on the left panel of Fig. 1111 
Finally, the relevant solution from a thermodynamical point of view is the one which maximizes V{m^ qi) over both m 
and qi ; because the continuous maximum is always very small and merges with the discontinuous one upon increasing 
m, one finds that the relevant solution is always the discontinuous one. 

Appendix B: 1-step quantum cavity equations and their numerical resolution 

We briefly recalled in Sec. Ill Al how the classical coloring problem on a diluted random graph could be solved 
analytically in the thermodynamic limit by the cavity method, summarized by the RS ^ and IRSB ^ cavity 
equations. The same route can be undertaken in the quantum case, with the important complication that the 
probability distributions rj appearing in these two equations are now over "trajectories" instead of colors, a trajectory 
(T being a piecewise constant periodic function from [0, f3] to {1, . . . , gcoi}- The only modification to Eqs. pi4p implied 
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by this change is the replacement of the functions g, z hy (see for more details) 




= rr"(-)^n\,(<T,)exp / 




(Bl) 



z{rii, . . . jf^c-i) being defined by normalization, and n^cr) counting the number of discontinuities in the trajectory cr. 
Note that for F = 0, only constant trajectories remain and one recovers the definition given in ([3]). 

The population dynamics method [sol [t^ is a convenient way to solve these equations numerically. Probability 
distributions are approximated by weighted samples of representative elements: 



We therefore have two deal with two level of populations: one of Afpop messages t], each of these messages being 
represented by A/traj trajectories cr. Each trajectory is encoded by the imaginary times in [0,/3] where it changes 
value (there are of order /SF such times), and by its constant value between these jumps. Then the IRSB equation 
(|4IB1I) can be solved by iteration on these samples and the associated weights , bij , the current estimation of P 
being inserted in the r.h.s. of Eq. Q, and the l.h.s. is approximated by a new discrete representation. This procedure 
would be exact only if Mpop and A/traj were infinite, while the memory available on present days computers limits 
the values of A/'pop and A/traj to rather small values (see below and App. [C] foi' concrete examples). This induces 
both systematic deviations of the empirical mean from the exact value and noise in its estimation; extrapolations to 
A/'pop, A/traj oo via finite size analysis can in principle be performed to reduce these effects, we show one example of 
such treatement in App. [CT] Further difficulties arise when the weights (or bij) become very heterogeneous: then 
the population is effectively supported only on the representants with the largest weights, hence the effective size of 
the population can be much smaller than the number of samples. These population representations of probability 
distributions are known in statistics as particle approximations (TSj : many resampling techniques are known to fight 
this impoverishment of the sample representativity, but there does not seem to be an universal way to avoid it. 

To increase the speed of our numerical code we made use of the parallelization opportunities offered by multi-core 
computers. Let us sketch how the IRSB equation resolution can be distributed on several processing units. The 
update procedure that we alluded to above can simply be thought of as a way to generate a sample at step t + 1, 
{{r]i{T + 1), ai(r + 1))}, from a sample at step r, {(?7i(T), ai(r))}. The important point is that the new representants 
{r]i{T+l), aiir + l)) are generated independently one from the other (apart from the normalization condition = 1 

that can easily be enforced). The update procedure can then be parallelized as follows: each core is first sent the whole 
sample of messages and weights {rii{T), ai(T)}i<i<j\f^^^ . Then each core c, with 1 < c < Aycoro, generates independently 
a new sample of A/'pop /Aycoro messages: {{r]c^i{T + l), ac,i(T + l))}i<i<A^p„p/A^^„„ • These new messages are then gathered 
together to form the full sample at step r+l: {(?yi(T), ai(r) )}i<i<7V'pop ■ This method does not change the memory hmit 
of the procedure, because each core is sent the whole population in the first step, but allows for a gain in time roughly 
proportional to A/"core (the communication between the processors usually takes much less time than the generation 
of the samples). Typically we used A/'core between 12 and 64; values of A/'pop and A/traj are limited by the amount 
of memory available on a single core, leading in our case to the constraint A/'pop x Atraj ^ 3.10''/(/3F). A moment 
of thought reveals that it is also possible to avoid the step of population gathering, and to keep the population of 
messages {rji^r), ai(T)}i<i<jv ^ spread over the A/'core cores. This allows for larger populations, but strongly increases 
the time needed to update the population, as much more information has to be exchanged between cores; therefore, 
we dit not use this second procedure in this work. 



The simplest transition to find numerically is the static continuous transition at T1(F). In fact, it corresponds to 
the point where the replica-symmetric maximum becomes a saddle point within the larger IRSB subspace; hence the 
RS solution, if slightly perturbed, will be unstable under iteration. Therefore, it is enough to initialize the IRSB 
population dynamics equation Q on the RS solution and to check whether this solution is stable under iteration, 
with respect to a small perturbation. This is detected by the growth of the intra-cluster overlap gi, defined in ^ 
in the quantum case, from its RS value I/qcoI- Formally expanding in (|4]) P around a delta-peaked distribution, one 
realizes that the condition of local stability of the RS solution towards a non-trivial IRSB solution is independent of 




(B2) 



Appendix C: Details on the numerical results and fitting procedures 



1. Finite population scaling for the instability of the RS solution 
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FIG. 12: Scaling form (ICT| for gT = gi - 1/gcoi, for c = 9, T = 0.09, and the small field instability (rr(r = 0.09) = 0.08) on 
left panel, large field instability {Tf{T = 0.09) — 0.44). The dotted lines are the linear asymptotics for the scaling functions 
F; have been adjusted to obtain the best data collapse. 



the Parisi replica-symmetry breaking parameter m. It is thus possible to take m = 0, which is very convenient in 
practice: the weights z{{rii\)™' are all equal to 1 (at non-zero temperature the factors z never vanish), the external 
population representation has thus homogeneous weights and A/'pop can be reduced to a rather small value without 
too much negative effects on the numerical accuracy. We can then concentrate on the finite population size effect of 
A/traj and build a scaling theory around the transition. To simplify the discussion we focus on the case in which the 
temperature T is fixed and T is lowered from the quantum paramagnetic phase at large T through the first continuous 
instability [T). First of all, the RS solution has a very simple overlap structure: because of the symmetry between 
colors, one has go = 9i = on the RS solution. Introducing qi = qi — l/(?coi, one finds that within the RS phase 
(r > rj^(T)), qi has a finite population size behaviour as 1/ ^/Mtvay On the other hand, we expect from the discussion 
of Appendix \K\ that for F < F^(T), qi — —a (F — Fj^(T)). These two behaviours can be matched by introducing a 
scaling function T such that: 

1 



gr(r,T,A/;,,j)--_^ (f -F+(r))/v;~ (ci) 

V-'^traj L J 

Moreover, T{x) admits the following asymptotic behaviours: lim2;_j._oo -^(a;)/a; = —a, lim^^-^oo -^(2;) = 0(1). The 
value of F — F;*" (T) is then determined in order to obtain the best collapse of the numerical data for different values 
of A/traj ; as shown in Fig. [T^] this scaling form gives a very good collapse of the numerical curves. This procedure 
therefore allows one to determine the line T'i(F) with very good precision. 



2. The low r regime for c = 9: the lines rK(r), T*(r), rd(r) 

Having discussed the determination of ri(F), let us now focus on the more complex determination of T'k(F) and 
Ti,(F) that were defined in Sec. IIII Bl and appeared in Fig. |4] and [5] We assume here that the reader is already familiar 
with the discussion of Sec. IIII BI but we will repeat some parts of the discussion for clarity. 

For simplicity, we will discuss the procedure at a fixed low temperature T = 0.09, and explain how Fk and F^, are 
determined at this temperature. At F = 0, there is a single IRSB solution of the classical cavity equation, the IRSBhq 
one, that exists in an interval m £ [rn^^, 1]. This solution can easily be found by solving the classical IRSB equations, 
it is the relevant one and has positive complexity at m = 1, because Td(F = 0) > 0.09 and TK(r = 0) = 0. The static 
Parisi parameter is therefore ms(F = 0) = 1. The evolution upon increasing F from F = has already been sketched 
in Sec. IIII Bl and in Fig. [H but we give some additional details here (see a summary on Fig. [T^: 

• For F < Fk(T), in addition to the RS solution, there exists only the IRSBhq solution, which is defined for 
m S [w^q(F), 1] and has a positive complexity at m = 1. It can be easily constructed by initializing the IRSB 
population dynamics in this solution at F = and following it at positive F. The value of "T-hql^) '^^^ ^® 
determined numerically by following the evolution of the IRSBhq solution upon decreasing m from 1 to at 
constant F, as the value at which the solution disappears and the complexity jumps to zero (see Fig. 1141 left 
panel). 
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FIG. 13: Domains of existence of the solutions of the IRSB equation for c = 9, T = 0.09, in the (F, m) plane. The IRSBhq 
solution exists for m > rn'^^{T), while the IRSBiq one exists for F > Fr and m < mjq(F). The (extrapolated) intersection of 
m^q(F) and mj^(F) defines the point F* at which the two IRSB solutions merge into a single one. The region of coexistence 
is therefore delimited by Fr < F < F* and rn^^{r) < m < mj^(F). However, only the IRSBhq solution is thermodynamically 
relevant: the value of ms(F) is equal to 1 for F < Fk, and it is smaller than 1 for F > Fk, but it always correspond to the 
IRSBhq branch. 



• For rK(T') < F < {T) there still exists only the IRSBhq solution, but its complexity is negative at m = 1. 
The value of T}<i{T) is therefore determined as the point where the complexity at m = 1 in the IRSBhq goes 
continuously to zero, see Fig. [131 The static value of the Parisi parameter ms{T) g ("^hq(^), 1) is determined as 
the point where the complexity vanishes (see again Fig. [Ml left panel). 

• For rj"(r) < r < r*(T) the RS solution becomes unstable and it disappears in favor of the IRSBiq solution 
which exists in an interval m G [0, "i-^q (F)]. Therefore in this region the two IRSB solutions coexist for m^q(r) < 
m < 77ijq(F). The IRSBhq solution can still be obtained numerically by following it in F at to = 1, and then 
decreasing m at fixed F. On the contrary, the IRSBiq solution is obtained by starting at to = in the RS 
solution; as the RS solution is unstable, the population evolves towards the IRSBiq fixed point which can then 
be followed by increasing m at fixed F. The point to^^ (F) is defined as the point where the population jumps 
into the IRSBhq solution. Again, in this region ms(F) £ (mjjq(F), 1) is the point where S = in the IRSBhq 
solution. 

• Upon increasing F, the coexistence region shrinks; the field F^(T) at which mj^q = TOjq marks the limit of 
coexistence of both solutions. However, when the coexistence region is too small it is hard to determine the 
spinodals, so the point F^, is obtained by extrapolating the two lines TOjljq(F) and TOjq(F) from smaller F to 
determine their intersection (see Fig. [T3|. 

• Finally above F* there is a single IRSB solution, and again TOs(F) is the point where its complexity vanishes. 

The evolution of the complexity in Fig. [Ml is the best way to provide numerical support to the schematic picture for 
the free energy (j){m) outlined in the right panel of Fig. [31 Unfortunately, a direct computation of (/)(m) is not possible 
because the data are too noisy. The fact that the relevant solution for thermodynamics is the IRSBhq one is therefore 
not evident from Fig. 1141 However, this can be deduced from the consistency arguments presented in Sec. IHIBI and 
from the fact that the complexity of the IRSBhq is much bigger, thus suggesting a more rapid rise of its free energy 
from the RS value. 

The dynamical transition temperature rd(F) can be determined quite easily as the point where the IRSBhq solution 
at TO = 1 disappears on increasing temperature at fixed F, like in the classical case. An example of this procedure is 
shown on Fig. [151 R is important to stress that when the two lines T'd(F) and rK(F) merge, the overlap qi does not 
tend to its RS value, as shown in Fig. 1161 This is an important difference with respect to the spherical p-spin model, 
where the two lines merge on a critical point where the transition becomes continuous. In the coloring problem, the 
overlap qi remains non trivial, indicating that the lines and Tk do not merge on the instability line Tj, consistently 
with the phase diagram reported in Fig. [5l 
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FIG. 14: Complexity E(m) for c = 9, T = 0.09 and T = 0.07 (left), T = 0.1 (middle), T = 0.2 (right). Results obtained 
with A/'pop = 1200 and Atraj = 3500. The red points are obtained by decreasing m from 1 down to 0, thus selecting the 
IRSBhq solution, while the black points are obtained increasing m from 0, thus selecting the IRSBhq solution (if any). Note 
the resemblance with Fig. [Til For this temperature, Fr ^ 0.064, Ff ~ 0.08, F* ~ 0.25. The left plot is for Fk < F < Ff : in 
this region there exists only the IRSBhq solution with E(m = 1) < 0. For Fr < F < F* (middle and right panel) there exist 
two branches IRSBhq and IRSBiq, which can be seen to merge when F approaches F*. To find the thermodynamic solution, 
one has to maximize the free energy of the system, which amounts to enforce the condition E(m) — (or E(m = 1) > 0): this 
condition always selects the IRSBhq solution. 
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FIG. 15: Complexity of the IRSB solution at m = 1, for the coloring problem with c = 9. Left panel: as a function of F for 
T = 0.09. Right panel: as a function of T for F — 0.13. No clear finite size effect could be found so we report the data obtained 
with three large sizes. The black line is a quadratic fit on the range F € [0.02, 0.12] (left panel), T € [0.14, 0.17] (right panel), 
giving Fk(T = 0.09) = 0.063 and Tk(F = 0.13) = 0.155. The value of Td is marked by an arrow. 



3. The spurious RS first order transition for c = 13: a study in the (F, m) plane 

As mentioned above, at the replica-symmetric level and for T < 0.12, the qcoi = 4, c = 13 model exhibits a first order 
phase transition when F is varied, while the transition is only third order when the IRSB computation at the static 
value ms(T, r) of the Parisi parameter is done. Altough the IRSB equation (U) yields thermodynamic observables 
only for m — ms{T, F), the static value of the Parisi parameter, one can also solve it for any value of m. This gives a 
way to interpolate between the first order phase transition at m = 1 and the continuous transition at m — ms(T, F). 
The results of cavity computations performed for various m are shown on the left panel of Fig. [T7]for T — 0.06. The 
value m — I corresponds to the replica-symmetric calculation, in which the hysteresis between the solution coming 
from F = and from large F is very well marked. Decreasing to, the hysteresis gets smaller and smaller, and cannot be 
seen anymore for to — 0.5. In this case TOs(0.06, ri(0.06)) ~ 0.13 and therefore the physical transition is continuous. 
As for the c = 9 model at small field, it is possible to build a scenario to explain these numerical results and to 
understand in greater details the nature of the transition; we briefly sketch it herafter. 

These results can be summarized by drawing a phase diagram in the (F, to) plane, for T fixed (right panel of Fig. lTTl) . 
At small F there exists only one IRSB solution, defined for all values of to. However for larger values of F there exists 
a range of m (between the blue and red lines of Fig. [T71 right panel) where two solutions of the IRSB equation coexist, 
similarly to what we explained in Sec. lIIIB II for q = 4, c = 9 model when F was decreased below F^,. In this case, one 
finds that the thermodynamic solution at mg always lies on the "low q" branch, defined for to £ [0, TOj^(r, F)]. Upon 
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FIG. 16: Difference between the intra- and extra-cluster overlaps on the Kauzmann transition line for the spherical 3-spin 
model (left panel) and the coloring problem with gcoi = 4 and c — 9 (right panel). The triple points he and To are defined in 
the text and do not have the same interpretation in the two cases. Note that the y-axis does not start from zero on the right 
panel. 
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FIG. 17; Coloring problem at gcoi = 4, c = 13 and T = 0.06. Left panel: Internal overlap qi{T, F, m) for several m as a function 
of F. The value of the continuous instability Ti{T) is indicated by an arrow. For m > 0.6 (including the RS case m = 1), 
hysteresis can be seen between the curves coming from high F and those coming from small F; but this hysteresis disapppears 
for m < 0.5, in particular for the static value ms{T,Ti{T)) ~ 0.13. Right panel: Phase diagram in the (F,m) plane. The 
replica-symmetry is broken on the left of the dashed black line, which represents the continuous instability of the RS solution. 
The high overlap solution exists above the red line with squares, and the low overlap one below the blue line with circles. These 
two lines meet at the value of m below which no hysteresis remain. The relevant solution is always the low overlap one; it 
becomes RSB for F below the continuous instability shown as a vertical dashed black line. Finally, we also report in this region 
the static value of the Parisi parameter ms{T, F). 



increasing the field, rn^ reaches 1, and above the value of the linear instabilty Ti{T), this low overlap branch becomes 
m-independent: this is the (quantum paramagnetic) replica-symmetric solution. On the other hand, the other "high 
q" branch still exists for F > Fj not too large, and m S [m^q(r, F), 1]. It completely disappears only when m^q(T, F) 
reaches I, which corresponds to the limit of existence of the corresponding replica-symmetric solution. This branch 
is never relevant, because maximizing over m leads to a solution of negative complexity, that has to be dismissed. 
However, it explains why a replica-symmetric computation gave the illusion of a first order phase transition. 



Appendix D: Perturbative expansions for small transverse field 

In this Appendix we explain how to derive a small transverse field perturbative expansion for observables in the 
quantum coloring problem. Let us first consider more generically an Hamiltonian H{T) — i/p + THq, where -ffp is 
diagonal in the computational basis, with diagonal elements denoted E{q_), and where Hq is purely off-diagonal. The 
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expansion in F of the partition function reads 



Z(/3,r)=Tr 
= Tr 















o(r3) 



1-5 



p{E{a!) - E{a)) 



0(r3) 



(Dl) 



This expansion can of course be obtained from the usual formulas for the perturbation of eigenvalues, taking into 
account their possible degeneracies. We now assume that i/q takes the form of an uniform "transverse field" — T^, 
with the flipping operators Ti defined in Eq. ([7|) . For simplicity of notation we shall work with the free energy density 
instead of the partition function. The expansion for this quantity reads 



N 



/(/3, r) = r = 0) - /?r2 ^ Ma)l ^ ^ (Is 



1-5, 



E{a}.E{g_li-'''>) 

1 a' ^G-i 



0(r3), (D2) 



where ct^*'*^ ^ denotes the configuration obtained from g_ by replacing its z-th variable by cr', and we introduced the 
classical Gibbs-Boltzmann probability law /^(ct) = e~^^^-'' jZiji^V — 0). Deriving with respect to F we obtain a 
similar equation for the transverse magnetization, 
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' l3{E{a(^''^')) - E{(7_)) 
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(D3) 



We would like to emphasize here that these first quantum corrections are expressed as an average of simple, local, 
observables with respect to the classical Gibbs-Boltzmann measure - unfortunately, for higher order terms in F, the 
corresponding classical quantities involve correlation functions between variables at arbitrary distances. In particular, 
if the classical energy E{a) is such that /i is correctly described by the classical cavity method (be it in its RS or 
RSB version), then the first quantum correction can be obtained by solving the classical cavity equations, which are 
much simpler than the quantum ones. As a matter of fact (|D2I) can be obtained from the quantum cavity method, 
either by expanding the quantum messages rj{a) in powers of the number of discontinuities they contain (0 for the 
classical trajectories, 2 at the lowest order in F) and noting that their first quantum correction can be expressed from 
the classical component of the incoming messages (see Eq. IB1|) or by taking the derivative with respect to F^ in the 
variational expressions of the free energy provided by the cavity method in terms of its order parameter (77 at the 
RS level, P{ri) at the IRSB one). In the IRSB case the quantum corrections are thus the sum over the postulated 
pure states of the correction inside each pure state; this tells indirectly something about the classical cavity method, 
namely that most pairs of configurations belonging to two different pure states are at Hamming distance strictly larger 
than 1, as it should. 

As usual in perturbation theory, the presence or absence of close-by configurations with degenerate energies leads 
to qualitatively different behaviors; this is particularly relevant at low temperatures, as the two corresponding terms 
in Eqs. (|D2ID3I) have different scalings with (3. In particular, for the coloring problem, the transverse magnetization 
given in Eq. (jD3[) satisfies 



lim lim = (gauth(2:,«) 
P^ooT^O /3F 



1) 



(D4) 



where in the r.h.s. the average is over the ground states g_ of the model, and qauth(o;, *) is the number of colors that i 
is authorized to take by its neighbors (in other words gcoi minus the number of distinct colors the neighbors of i take 
in a): the more "floppy" the spins are, the stronger they respond to the quantum field. As mentioned above these 
small F expansions can be obtained from the quantum cavity formalism, at any level of replica-symmetry breaking. In 
particular at the IRSB level the formula (jD4p can be interpreted as an expansion for the transverse magnetization of 
a single cluster, provided the average in the r.h.s. is restricted to this particular classical pure state. This provides the 
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justification for the interpretation of the Quantum Monte Carlo simulations presented in Sec. IIV Al and in particular 
of Eq. (HH). 
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